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ABSTRACT 


The two-dimensional (2D) spectral properties, including the degree of isotropy, of 
under-ice roughness in the Beaufort Gyre were investigated. Under-ice roughness plays 
an important role in turbulent ocean heat, salt, and momentum fluxes that determine the 
delicate balance between surface forcing and the ocean interior that sustains or reduces 
the perennial ice cover. Accurate characterization of this roughness is important for 
numerical modeling and prediction of the Arctic air-ice-ocean system, which will play a 
significant role as the US Navy increases its strategic presence in the ice-diminished 
Arctic Ocean. Ice draft data from moored upward looking echosounders of the Beaufort 
Gyre Exploration Project and ice motion derived from trajectories of International Arctic 
Buoy Program buoys were combined to create directional sections of ice draft. 
Autocorrelation functions were created for each segment which were then combined to 
make a 2D autocorrelation function. The 2D FFT was taken of this in order to yield a 2D 
spectrum, from which isotropy and other spectral properties could be determined. The 
processing routine was tested on several synthetic test datasets and in each case produced 
the correct output. Two years of data from three moorings in the Beaufort all provided 
2D spectra which showed weak anisotropy of the under-ice roughness. Error estimates 
were determined by Monte Carlo simulations. For each dataset, most of the variance of 
the underside morphology was contained in a wavenumber band between 0.001 and 
0.01m" (100-1000m). There were no distinct, preferred wavelengths. Spectral levels 
ranged from 1.2 - 4.2 x 10* m’ for Mooring C (2004-2005) to 0.2 — 1.6 x 10° m’ for 
Mooring A (2004-2005). 
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I. INTRODUCTION 


A. CURRENT EVENTS 


Recent changes in the Arctic Ocean and specifically its diminishing ice cover 
have received significant attention both in scientific circles and amongst the general 
public. The current International Polar Year is bringing increased awareness and 
gathering crucial data needed to understand these changes. An ice-diminished Arctic is 
also of strategic significance to the national security of the United States. The viability of 
both the Northern Sea Route and the Northwest Passage as navigable shipping lanes for 
the first time in recorded history combined with the new feasibility of resource 
exploitation (minerals, oil, fisheries) in the Arctic Basin itself has led many nations to 
reassert or even stake new claims in the region. Russia, Norway, Sweden, Finland, 
Denmark, Iceland, Canada and the United States all have major interests in the Arctic and 
each seeks to extend their Exclusive Economic Zone (EEZ) in any way possible. A 130 
year long dispute continues between the United States and Russia over the location of a 
line dividing possession of the Bering Sea between Siberia and Alaska (Damron, 1997). 
Norway and Russia have competing claims to fishing rights in the Barents Sea (BBC 
News, 2005). Canada insists that the Northwest Passage lies within Canadian territorial 
waters and that they will assert sovereign control of it whereas the United States and 
many other countries claim international rights to use it for shipping (Shukman, 2007). 
Most recently, in the summer of 2007, the Russian Federation submitted a formal claim 


for over 45% of the Arctic seafloor, including the North Pole itself. 


The summer of 2007 also saw several military developments. A Russian 
submersible symbolically planted a flag on the ocean floor at the geographic North Pole. 
Canada announced that it would enhance its military presence in the Arctic, including 
building two new military facilities; one for cold-weather warfare training at Resolute 
Bay and the other a deep-water port on the northern end of Baffin Island (Reid, 2007). 
The Canadian government also placed an order for six to eight new Polar Class 5 Arctic 
Offshore Patrol Ships. The vessels will be armed and according to Prime Minister 


1 


Harper, they will “allow us patrol the Arctic waters in the Northwest Passage when it 
matters (CBC News, 2007).” An American guided missile cruiser, USS NORMANDY 
(CG 60), conducted a patrol across the Arctic Circle north of Iceland to show that the US 
Navy maintains the capability to operate in the far north. USCGC HEALY conducted a 
bathymetric survey of the Chukchi Cap in coordination with NOAA’s Office of Coast 
Survey “to better understand its morphology and the potential for including this area 
within the United States’ extended continental shelf under the United Nations Convention 
on Law of the Sea (Elliott, 2007).” In this “gold rush” atmosphere, there is developing 
potential for armed conflict over both resources and access. At the very least, the US 
Navy can be expected to conduct routine Freedom of Navigation demonstrations, sending 


ships to physically challenge territorial claims that the US does not recognize. 


In order to safely operate in the Arctic, whether a vessel is military, merchant, or a 
tourist cruise ship, ice modeling and the resulting forecasting will become as important as 
routine weather forecasting is today. There are a number of challenges in ice modeling. 
Due to the remote and harsh Arctic environment, these difficulties include sparseness of 
in situ data. Furthermore, there is a limited amount of remotely sensed data. Of all of the 
unknowns, one that is particularly challenging is morphology of the underside of the 
icepack which actually is in contact with the ocean. This cannot be determined remotely 


with current technology. 
B. UNDER ICE ROUGHNESS FEATURES AND DEFINITIONS 


Pressure ridges are formed in the Arctic Ocean when two ice floes are forced 
together by wind or current. During the collision, areas of thin ice between the floes are 
crushed and sometimes even the edges of the floes themselves fracture. This material 
accumulates and piles up so as to form a ridge visible from above the ice. An 
approximate mirror image (but much larger in size) also forms as an inverted ridge 
protruding downwards into the ocean away from the underside of the ice floe. This 
submerged feature is referred to as an ice keel. It has been estimated that pressure ridges 
and their associated ice keels account for approximately half of the ice volume in the 


Arctic Basin (Davis and Wadhams, 1995). 


C. SIGNIFICANCE 


Ridging events influence the magnitude and direction of ice motion. Pressure 
ridges cause an increase in the aerodynamic drag at the air-ice interface and their 
associated keels cause an increase in hydrodynamic drag at the ice-water interface. If an 
ice keel is deep enough it can generate internal waves at the pycnocline, exerting much 
higher drag forces. This can lead to the phenomenon referred to by Ekman as “dead 
water” (Wadhams, 2002). This occurs in highly stratified areas where much of the 
energy used by a ship for propulsion is lost due to the formation of internal waves as the 
keel plows though the pycnocline. This loss of momentum to internal wave generation 
complicates modeling of the total air-ice-ocean momentum budget, especially when the 


size and fregency of the keels is unknown. 


It is also important to note the flow effects induced by ice keels. The ratio of a 
keel depth to the depth of the oceanic boundary layer is very large and can even exceed 
one, as was the case at the Applied Physics Laboratory Ice Station (APLIS)-2007. This 
makes a field of large ice keels analagous to a large mountain range in the atmospheric 
boundary layer. Keels can cause flow acceleration, blocking and stagnation along with 
turbulent fluxes of heat, salt, and momenturm. Simulations indicate that the turbulent 
wakes generated as ice keels are pushed through the water can exend several hundred 
meters downstream. Large Eddy Simulation modeling has shown that a 10m deep keel 


can cause a five-fold increase in heat flux (Skillingstad et al, 2003). 


There are engineering applications which need to account for ridging. Ice ridges 
are the deciding factor in determining strength requirements for Arctic offshore platforms 
and icebreaker hull design. The depth of ice keels and their potential to scour the bottom 


in shallow water determines the depth of pipeline burial (Wadhams, 2002). 


There are significant naval applications as well. Safety of submarine navigation 
in the shallow waters of the Arctic shelf and ability of a submarine to surface through the 
ice are affected by ice keels. Possibly even more important are the acoustic effects 
attributed to ice keels. Due to its density structure, the Arctic Ocean is characterized as 


an upward refracting acoustic environment (Fig. 1). Acoustic energy gets turned towards 
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the ice where it hits and is scattered as a function of the ice bottom roughness (Fig. 2). 
An even more extreme case is where an ice keel is deep enough to cause blocking of the 
acoustic path (Bourke and McLaren, 1992). The acoustic shadows cast by an extremely 
deep keel make it possible for a submarine to hide in its lee, shielding it from active 
sonar. Also of significance to submarine operations is the ambient noise of the ice pack 
itself. High levels of ambient noise are associated with ice floe collisions, especially at 
the ice edge (Wadhams, 2002). This phenomenon could greatly complicate passive 


sonar detection. 
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Figure 1. Simplified Arctic sound speed profiles. Model | is based on a sequence of 
plane parallel layers each with a constant gradient whereas Model 2 is based on a series 
of flat layers of constant velocity and density (Kutschale, 1969). 
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Figure 2. Acoustic ray paths for a source at 100m depth using Model | from Fig. 1 
(Kutschale, 1969). 


The effects of ridging may extend to mechanisms important to climate change. 
The more ridged the ice is, the greater the volume of fresh water that can cross flux gates 
such as Fram Strait. Ridging also increases the melting rate of ice, by increasing the 
surface area of the ice in contact with the warmer ocean underneath. If an ice keel is 
deep enough, it can cause deeper (warmer) water to be dragged up behind it and 
subsequent turbulent mixing will bring this warm water into contact with the underside of 
the ice, further melting the ice from below. Observational data support the assertion that 


ridged ice melts faster than undeformed ice (Davis and Wadhams, 1995). 
D. STUDY OBJECTIVES 


The purpose of this study was to produce a two-dimensional (2D) spectral 
respresentation of ridge properties and roughess of the underside of Arctic sea ice which 
would allow determination of anisotropy. The spectral approach makes it possible to 
determine preferred direction, wavelength, and magnitude of the roughness. If the 
roughness is isotropic, coupled air-ice-ocean models can reliably use a single coefficient 
for roughness, regardless of the direction of the wind-forced ice motion. However, if the 
roughness is anisotropic, directionally variable, direction-dependent coefficients may be 
necessary to accurately simulate ocean-ice exchange rates. This study combines moored 
upward looking sonar (ULS) derived ice draft data with motion determined from drifting 


buoys to sparsely sample and develop 2D spectra of under-ice roughness topography. 
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HW. DATA SELECTION 


A. ICE DRAFT DETERMINATION 


Significant efforts have been made to infer under-ice topography from airborne 
laser techniques. There have been simultaneous submarine upward looking sonar and 
and LIDAR equiped aircraft transits that attempted to establish relationships between the 
two surfaces. While there was some success, the amount of snow cover is variable and 
difficult to determine, thus contaminating the freeboard to draft ratio for the ice alone 
(Wadhams, 2002). A typical ratio used as a “rule of thumb” for ice keel draft to ridge 
height is 3.2:1. Recently, measurements of a ridge at APLIS 2007 did not follow the 
relationship. A hole was melted through the ice near an ice ridge to deploy instruments. 
The ice was approximately 13.5m thick, which suggests that the hole went through the 
ice keel itself. The associated pressure ridge had a maximum surveyed height of 2.88m 
above sea level, implying that an its underlying keel would be on the order of 9.2m depth 
(Fig. 3). Even if the tip of the keel had been offset from the ridge through some sort of 
shearing or rafting event, it still would not closely follow the ratio. Furthermore, any 
mirror image of ridge and keel geometry would to some extent be distorted for multiyear 
ice, which has undergone annual cycles of melting and freezing cycles, along with 
mechanical stress and strain events and their subsequent deformation. This observation 


shows significant differences between ice keel and surface ridge structure. 
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Figure 3. Surveyed ice ridge from APLIS 2007 and an extrapolated ice keel using 
the stadard ratio for ice keel draft to ridge height of 3.2:1. 


There are two common methods of obtaining reliable ice draft distribution data, 
both of which use ULS. One method is to mount an ULS on a submarine which then 
operates under the ice, recording ice draft as it goes. The second method is to mount the 
ULS to a bottom mooring which extends to within about 50m of the surface and record 
the ice draft which comes into its field of view. There are strengths and weaknesses in 


both methods. 
1. Submarine Data 


The strength of submarine data is that it can give a near-instantaneous snapshot of 
a broad area because the submarine can move quickly in relation to the speed of ice 
motion. During submarine missions, ULS data storage and power requirements are not 


major concerns. 


However, submarine data have a number of limitations. Data can be 
contaminated by vessel depth or speed variations. Submarine cruises are spatially and 
temporally sparse. The areas surveyed are rarely if ever repeated. Furthermore, except 
for the SCICEX missions, collecting environmental data is not the submarine’s primary 


mission and therefore quality may suffer somewhat. 


Submarine ULS data are rarely released publicly due to security issues concerning 
sovereignty of territorial waters. Submarine datasets were rarely released to the public 
before Al Gore, then Vice President, petitioned to declassify data in an area which has 
become known as the “Gore Box” (National Science Foundation, 1998). This release 
area has since been expanded by the Chief of Naval Operations beyond the original 


“Gore Box” (Fig. 4) (National Snow and Ice Data Center, 2006). 





Figure 4. Original “Gore Box” and current release data for some submarine 
collected upward looking sonar ice draft data (National Snow and Ice Data Center, 2006). 


There are additional safety issues to consider. It is dangerous for a submarine to 


be in shallow water under the ice with little room to maneuver. For example, USS 
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SARGO (SSN 583) damaged her sail in a collision with an ice keel in the shallow and 
confined waters of the Bering Strait during the winter of 1960 (Fig. 5) (Nicholson, 2008). 





Figure 5. Crewmembers inspect sail damage on USS SARGO (SSN-583) after 
colliding with an ice keel in the Bering Sea in February 1960 (Nicholson, 2008). 


Zz Mooring Data 


Mooring mounted ULS allows for the collection of a long term record of sea ice 
draft. The locations are usually re-occupied so multiple years of data can be collected in 
this manner. It is cheaper and safer to deploy a mooring than to operate a nuclear 
submarine. Moorings can be located in the very shallow waters of the continental shelf 


where submarines cannot operate. Finally, these types of datasets are not classified. 


Moored ULS has weaknesses too. Temporal coverage is good but spatial 
coverage is limited to a single point for each mooring. The unit must be completely self- 
contained and rugged enough to survive the harsh conditions, placing restrictions on 
memory capacity and power requirements. Another weakness is that the data cannot be 
examined during the period of collection. The failure of a unit will not be apparent until 


recovery, which is typically a year after deployment. 
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After comparing the available submarine-mounted and moored ULS data sets, it 
was decided to use the moored ULS data because the observations were recent, year- 
round, and available for a number of different locations. These attributes make seasonal, 
annual, and regional comparisons of under-ice roughness possible. Specifically, the data 
used in this project were from the Woods Hole Oceanographic Institute’s Beaufort Gyre 


Exploration Project (BGEP) (Proshutinsky, 2008). 
B. CONVERTING TO A SPATIAL SERIES OF ICE DRAFT 


The BGEP moored ULS provide time series of ice draft, z(t). However, 
determination of isotropy requires a 2D spatial determination of ice draft, z(x,y). To 
convert from one to the other, the time series is multiplied by the sampling rate and 


advection speed of the ice. Therefore, the motion of the ice needed to be determined. 
C; ICE MOTION DETERMINATION 


Ice motion can be derived from a variety of methods. These include Eulerian 
measurements with an acoustic Doppler current profiler (ADCP) pointed upwards from a 
fixed location, Lagrangian motion determined from drifting buoys, and satellite remote 


sensing products. 


Ideally, there would be an ADCP collocated with every ULS to give high 
resolution local ice motion. Unfortunately, this type of data is not always available. For 
this project, ice motion was determined from the position fixes of buoys coordinated by 
the International Arctic Buoy Program (IABP) (Rigor, 2002). This buoy-derived motion 
was preferred over a number of satellite products such as those derived from AVHRR, 
scatterometry, and other merged datasets because of its higher temporal and spatial 
resolution (Fowler, 2003). Most satellite products give averaged 12 or 24 hour motion 
fields. However, in order to measure 2D under ice roughness, the motion vectors had to 
be available with a temporal resolution on the order of one hour so that the frequently 
dominant ice motion at these latitudes, inertial oscillations, could be resolved. This high 


resolution data could only be obtained from the IABP buoys. 
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Il. DATA PROCESSING 


A. ICE DRAFT DATA PROCESSING 


As part of BGEP, three ULS units were mounted on fixed moorings (referred to 
as ‘A’, ‘B’, and ‘C’) under the ice of the Beaufort Gyre (Fig. 6). 


Mooring Locations for BGEP 


2003 - 2005 
mn wae, 





Elevation (m) 





Dist. from N. Pole on 180° Lon. (m) 





* ts) $ “10 | -15 -20 
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Figure 6. BGEP Mooring Locations. 


The moorings were anchored in water over 3500 m deep, with the uppermost 
instrument on the mooring being the ULS, which was located approximately 50m below 
the ice bottom (Fig. 7). The ULS instruments used were ASL Environmental Sciences 
model IPS-4, which operated at a frequency of 420 kHz with a beamwidth of 1.8 degrees, 
yielding a footprint on the underside of the ice of about 2m. Range between the ULS and 
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the ice bottom was sampled at 2 seconds. Total hydrostatic pressure and water 


temperature were sampled every 40 seconds (Krishfield and Proshutinsky, 2006). 
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Figure 7. Schematic of ULS sees (Krishfield and Proshutinsky, 2006). 


In general, ULS data require significant processing and despiking before they are 
usable. For this project in particular, ice draft was calculated by correcting the measured 
range for tilt and sound speed variations and then correcting the total hydrostatic pressure 
measured at the transducer for variations in atmospheric pressure and changes in seawater 
density. Sound speed and density are only determined from CTD casts at the beginning 
and end of the deployment. During this deployment, it was impractical to have a CTD 
continuously profiling so other methods are required. Temperature measured at the ULS 
can locally correct sound speed and density, while corrections can also be determined by 
using the range observations during periods of open water (where draft=0). However, 
this method has limitations because it is often difficult to distinguish frazil and thin ice 
from open water. Since in situ observations of atmospheric pressure are very sparse in 
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the polar regions, they were obtained from NCEP sea level pressure reanalysis data. 
After all corrections have been made, draft calculations have an error of +5—10cm. 


(Krishfield and Proshutinsky, 2006). 


A large source of error in all ULS observations is the finite beam width of the 
transducer. Although it is assumed that the draft value assigned to a particular location 
on the ice is along the transducer’s axis, this is not always the case. On sloping features, 
the first return of the sonar pulse may be due to an off axis feature which is shallower. 
This will cause the apparent draft at a point to be greater (Fig. 8). The net effect is to 
overestimate ice draft and to broaden the apparent width of features (Melling et al, 1995). 
The only way to mitigate this effect is to use an ULS with as narrow a beam width as 
possible. The IPS-4 transducer has a beamwidth of 1.8 degrees which is narrow in 
comparison to sonars used to generate previous datasets, which used beamwidths on the 


order of 14 degrees (Wadhams, 1978). 


APPARENT PROFILE 


TRACK OF 
SUBMARINE 





Figure 8. Effect of echosounder beam width on apparent keel profile. Note that the 
apparent width and average draft increase (Wadhams, 1978). 


Other sources of error are due to phenomena occurring in leads. Waves can give 
an apparent draft if sampling occurs while its trough is over the ULS. Likewise, the crest 


of a wave can give a solution with a negative draft, which would physically mean the 
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sonar was looking out of the water. Another effect is from wave breaking, which can 
inject clouds of bubbles into the water. These in turn may appear as legitimate ice drafts 


in the sonar record (WCRP, 2004). 


Another major complication in interpreting ULS observations is the effect of 
pressure ridge orientation relative to the direction of sampling. If a keel is sampled 
perpendicular to its axis of lineation, there is no problem. However, if the keel is 
sampled at some other angle, the effect is that the keel appears to be broadened and the 
spacing between features of similar alignment is lengthened as well (Fig. 9 and 10). The 
wavelength will always appear lengthened; it cannot appear shortened. This leads to a 
“reddening” of the spectra whereby actual wavenumbers are shifted to lower 
wavenumbers. The most extreme case occurs when sampling along the axis of a keel, 
where the wavelength or spacing between keels becomes infinite and wave number goes 


to zero. 
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Ice Keel with Two Sampling Directions 
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Figure 9. (a) View from below ice keel showing sampling directions perpendicular 


to the keel and cutting across it at 75°. (b) Apparent keel profiles as a function of 
sampling angle. @=0° and 6=90’ correspond to sampling directions perpendicular and 
parallel to the keel axis, respectively. Note that the apparent width and average draft 
increase as the angle becomes greater. 


17 


Field of Ice Keels with Two Sampling Directions 
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Figure 10. (a) View from below field of keels with sampling directions 0 and 75°. (b) 
Apparent draft profiles from field of regularly spaced keels as a function of sampling 
angle. Note the increase in apparent spacing of the ice keels as the angle increases. 
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This problem, sometimes refered to as spatial aliasing, needs to be resolved in 
order to generate a meaningful 2D spectrum of under-ice roughness. As the examples 
(Fig. 9 and 10) above have illustrated, a 1D slice through a 2D field is reddened (except 
for the one case where the slice happens to be exactly perpendicular to the orientation of 
the field of features). This reddening occurs because the 2D wavenumber information is 
lost. If the whole 2D spatial field were available, obtaining a correct 2D spectrum would 
be easy with a 2D FFT. This project used a processing routine for sparse data in which 
directional information was retained in the autocorrelation function, which preserved the 


2D wavenumber in the spectrum. 


B. ICE MOTION DATA PROCESSING 


Numerous IABP buoys passed through the area of the BGEP moorings from 
2003-2005 (Fig. 11). The associated raw buoy position data were obtained from the 
IABP office. Each file represented one buoy and covered a period of about | year. The 
following processing routine was applied to each data file using an algorithm written in 


MATLAB. 
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Drift Tracks for All Buoys 
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E 

° 1000 
e os 
3 
ad 
3 ry 
a _ 
c | E 
o 1000 $ 
= : 
° > 
@ 15 a 
5° 2000 
° 
= 2 
6 3000 
= 





Dist. from N. Pole on 90° Lon. (m) , 45 


Figure 11. Drift tracks of all buoys in mooring area from 2003-2005. 


First, the IABP buoy data file was loaded and all times were converted to year day 
referenced to 01 January 2003 as the start point. Next, redundant records were removed. 
The record was then scanned to find times when neither the latitude nor longitude 
changed for successive times. These records were removed as well. Position data were 
then converted from geographic coordinates of latitude and longitude to a polar 
stereographic grid compatible with the International Bathymetric Chart of the Arctic 
Ocean (IBCAO). This coordinate system is centered with the origin at the geographic 
North Pole. The x-axis measures distance in meters from the North Pole along the line of 
90 degrees East longitude and the y-axis measures distance in meters from the North Pole 
along the line of 180 degrees longitude. The buoys did not report position data on a 
standard time schedule. Therefore, the data were interpolated onto a common hourly 
time vector from the start of 2003 to the end of 2005. This allowed buoy positions to be 
available for every hour (on the hour). The polar stereographic buoy positions were 
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despiked to remove anomalous position data. A remaining problem was occasional, long 
temporal gaps. Interpolating across these would have introduced erroneous velocity 
estimates. These gaps in the position time series could be due to rime building up on a 
buoy’s GPS or broadcast antenna at some point and then getting washed off later. 
Significant data gaps were identified by finding position fixes that occurred more than 6 
hours apart. Interpolated data between these fixes were discarded. Velocity components 
7% 


(u,v) were then caluculated as a finite difference using u =—= 


The speed 
At tt 


(magnitude of the overall velocity) was determined via the Pythagorean theorem. 


The next task was to use the processed buoy data to generate a time series of u 
(east-west velocities) and v (north-south velocities) at each mooring. The following 


process was carried out for each mooring location. 


First, the mooring’s geographic coordinates of latitude and longitude were 
converted to the polar stereographic grid previously described. For each time step (every 
hour), the distance between all contemporaneous buoy positions and the mooring were 
calculated. All buoys within 250km of the mooring (Fig. 12) at that time were flagged 
and the mean of their velocity components (u,v) was then recorded as a value 
representative for that mooring. A time series was created at each mooring by repeating 
the process: advancing one time step at a time, determining the buoys within 250km, and 


recording the mean of their motion data. 
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Mooring C 
* All buoys within 250 km 
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Figure 12. Tracks of all IABP buoys passing within 250km of mooring C from 2003- 
2005. 
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Ice Motion Time Series 
Mooring C: August 2004-August2005 
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Figure 13. Example time series of motion compiled at Mooring C: August 2004- 


August 2005.need to correct 


Note that the actual direction of ice motion for mooring C was constantly 
changing because of the inertial oscillations (Fig. 13). Therefore, the motion data needed 
to be segmented into a number of segments so that each piece was approximately 
unidirectional. For the entire time series of motion (at each mooring), the direction was 


; ; afv 
calculated from u and v using a four quadrant arctangent function @ = tan ' (=| . Hourly 
u 


change in direction as well as cumulative change in direction were derived. Also, the 
hourly distance travelled and the cumulative distance travelled were calculated. If the 
cumulative change in direction or the cumulative distance travelled exceeded a threshold 
(defined here as 45 degrees for direction and 2,000m for distance), the segment would 
end with the point previous in time, so as not to exceed the threshold. The cumulative 
counters were reset to zero and the process was repeated. Segment start and stop times, 


and average speed and direction were recorded. 
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A limitation of this approach is that it assumes that the ice motion is continuous in 
space between the buoys and moorings, ignoring spatial gradients of velocity. For 
combining mooring and motion data, Melling et al (1995) cautioned against using 
Lagrangian methods such as radio beacons or GPS receivers in buoys mounted in the ice 
because they quickly leave the moorings behind and do not give true motion at the 
mooring. The error is particularly large near coastal boundaries and during the ice melt 
period. In this case, well away from land boundaries, it is assumed that the primary 
motion is due to inertial oscillations resulting from synoptic scale meteorological events 
of 1,000-2,500 km. With forcing of this scale, it is assumed that ice motion is continuous 


for a distance of 250 km. 
C. COMBINING THE ICE DRAFT AND MOTION DATA 


Without motion data, all that is recorded at the ULS is a time series of ice draft. 
The addition of motion data helps to transform this time series into a 2D spatially 
referenced description of ice draft. The process used below was first outlined by Hibler 


(1972). 


For each mooring, the ice draft information was loaded and corresponding times 
were converted to year days referenced to 01 January 2003 to be consistent with previous 
data. Using each time segment derived from the motion data, the corresponding ice draft 


time series had its mean removed and the autocorrelation function was calculated using 


N-1l 
R, (Ax) = aoe f (3) f (x, + Ax), where R,, is the autocorrelation value (in units of m ), 
i=0 


N is the total number of data points, f(x,) is the value of the draft function (in units of 


in units of m) evaluated at x (in units of m), and Av is the spatial lag (in units of m) (Fig. 


14). 
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Figure 14. Ice draft time series, distance series, and autocorrelation function for 
Segment #50 at Mooring C. 


The spatial lags Ax were a function of segment velocity. These lags were 
calculated using Ax=At-Vu>+v~, whereAt was the sampling interval (2s), and 
Vu’ +v’ was the segment average speed (m/s). 


For each mooring, all segment derived autocorrelation functions were then plotted 
together as functions of lag and segment direction. The result was a series of 
autocorrelation functions emanating radially from the mooring like spokes of a wheel 


(Fig. 15). 
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2-D Autocorrelation Functions 
Mooring C: Aug 2004-Aug 2005 
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Figure 15. Combined autocorrelation spatial coverage. 


It is important to note that the process used here departed from the one outlined 
by Hibler. Hibler recommended normalizing the autocorrelation functions. However, 
normalization strips away the total variance, which is critical for interpreting spectra. 
Normalized autocorrelation functions are unable to tell the difference between small and 
large roughness elements. The large elements, the ice keels, need to be distinguished. 
Their effects are distinct from roughness elements of a much smaller centimeter scale. A 
field of small roughness elements oriented in one direction could dominate a normalized 
2D spectrum, even though there was a more sparse field of very large roughness elements 
at a different orientation which in reality dominates the physical processes. Hence, in this 
study, variance was conserved (i.e. the autocorrelations were not normalized) in order to 


make direct comparisons between different locations, times, and directions. 


It is also important to note that this step of combining the autocorrelation 
functions is crucial because it mitigates the effects of sampling across features at an 
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angle. The relative change in autocorrelation functions between high and low lag 
numbers is determined by the direction that the linear ridge feature is crossed. By using 
the autocorrelation function and referencing it with its direction, the subsequent 2D 
spectrum does not become reddened because the 2D wavenumber is conserved. Energy 
peaks will occur at the correct wavenumbers even though the field is sparsely sampled 


This is demonstrated with test cases in section III A. 


This combined autocorrelation function is an intermediate step. The desired final 
product is the 2D spectrum. Since the autocorrelation function and the spectrum are a 
Fourier transform pair, it is straightforward to go from one to the other using the Fast 
Fourier Transform (FFT). However, the data needs to be on a continuous grid in order to 
use the FFT. The autocorrelation function had some gaps in coverage and so needed to 


be processed to be ready for the 2D FFT. 


The next step was to combine all of the autocorrelation functions by binning 
them in polar coordinates. The bins chosen were 10 degrees wide and 5m long to achieve 
satisfactory point density. It became apparent that each autocorrelation function had to 
be interpolated onto a common vector in order to keep the spacing of the points constant 
when compared to the other functions. It was observed that due to the fixed sonar 
sampling rate of 2s, autocorrelation functions corresponding to periods of slow ice 
motion had their data points spaced much more closely together than did functions from 
periods of faster ice motion. The result was that a particular bin might have 10 points 
from the slower autocorrelation function and only 1 point from the faster one, even 
though they both crossed the bin at the same angle. The value subsequently assigned to 
the bin was biased in favor of the slower moving segment of data. The autocorrelation 
data was therefore interpolated onto vectors of common length and point spacing, 
specifically 2000m and 0.5m respectively. This issue of oversampling during periods of 
slow ice motion and undersampling during periods of fast ice motion was addressed by 
Melling and Riedel (1995) who showed that resampling at a constant spatial rate was 
more appropriate than using the time series when trying to derive independent data for 
statistical comparisons. The mean was taken for each bin and its value was assigned to a 
point at the bin’s center. These binned data were then interpolated and transferred from 
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the polar coordinate system to a Cartesian grid measuring from -1300 to 1300 m in both x 
and y directions with grid points every 5m. The 2D FFT of the 2D autocorrelation 


function was taken to obtain a 2D spectra using 
1G 2 i(k, Ax, +1, Ay ;) : F 
S(k,1) = He aeRO, Ay pe i)’, where k and / are the horizontal and vertical 
i=0 j=0 
wavenumbers, Ax, and Ay, are the lags in the x and y directions, and R is the value of the 
binned and interpolated 2D autocorrelation function described above. 


The highest resolvable wavenumber, the Nyquist wavenumber was 


Nyquist ae, where Ax=grid spacing (5m). This corresponded to a wavelength of 


10m, which was expected to be suffucient to characterize the rapidly changing portion of 


the under ice environment. 


Synthetic data was created from a known field in order to test all phases of the 
process (see section III for detailed results). The process showed itself much more 
sensitive to variations in angular rather than radial bin width. Decreasing angular 
resolution caused the spectra to get smeared out. Increasing it caused critical decreases in 
point density in the inital binning process due to the sparse nature of the data. Radial 
variations had little effect on the spectra. However, as the radial resolution increased, 
there was not enough point density to fill every bin. If it decreased, the ultimate 
resolution of features occuring with higher wave numbers was decreased and roughness 
on these scales cannot be determined. In the end, both angular and radial resolutions 
used were a compromise in the fact that they had to remain coarse enough to deal with 


the sparse data. 


The size of the 2D autocorrelation area taken for the 2D FFT was determined by 
the length scale of the phenomenon occuring with the lowest wave number. This 
phenomenon was assumed to be distinct ice keels. Historical data shows the spacing 
between distinct ice keels has over a 99% probability of being less than 1000m (Fig. 16) 
(Wadhams and Davy, 1986). Therefore, the output from the FFT would need to be able to 


resolve a wave number of 0.001m™. Using resolution = = it follows that for a field 
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with grid spacing Ax equal to 5m, the required resolution of 0.001m' can only be 
achieved if the number of points N is greater than or equal to 200. This forces the size fo 
the field to go into the 2D FFT to have a radius of at least 1000m and the field to measure 


at least 2000m across. 
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Figure 16. Cumulative distribution functions of ice keel spacing. Plot (a) is derived 
from narrow beam sonar in the Beaufort Sea: open circles have a threshold of Om and 
closed circles have a threshold of 3m. Plot (b) is derived from wide beam sonar in the 
Eurasian Basin: open circles have a threshold of 0m and closed circles have a threshold 
of 15m (after Wadhams and Davy, 1986). 
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Even when using the synthetic data of idealized cases, there were artifacts in the 
spectra from the processing routine. To attempt to decrease these artifacts, zero mean 
random noise was added to the 2D autocorrelation surface before the 2D FFT converted 
the field into spectra. The noise level used was 1% of the standard deviation of the entire 
2D autocorrelation surface. Even with this improvement, there continued to be noise at 
very low wave numbers corresponding to long wavelengths. 2D convolution filtering 
failed to make the plots easier to interpret. Much of the noise was at wave numbers 
below 0.001m", which had previously been set as the maximum length scale imposed by 
independent keel spacing. Therefore, although it was apparent that the process could not 
differentiate between the processing artifacts and real data at very low wave numbers 
(long wavelengths), it was considered to be acceptable because energy is not expected to 


be there due to the physical scaling argument. 
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IV. ANALYSIS AND DISCUSSION 


A. TEST CASES 


To confirm that the technique worked as expected and to aid the interpretation of 
real world data, a number of test cases were conducted where a known field of under ice 
features was put through the processing routine. The random fields were generated and 
then sampled using real ice motion observations. The segments used were from Mooring 
C from August 2004-August 2005. For consistency, the Mooring C trajectories were 


used to generate ice draft data from the known fields for all of the test cases. 
1. Single Sine Wave 


The first test case was that of a field created from a single sine wave with 
amplitude 5m, wavelength 200m, and oriented so the crests form a line running vertically 
across the field (Fig. 17(a)). When put through the processing routine, the 
autocorrelation function showed the orientation correctly (Fig. 18), and the resulting 
spectra put the peak energy at the correct wave number and direction, [k,/]= [0 


(Fig. 17(b) and (c)). 
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Field from Single Sinusoid 
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Figure 17. (a) Field generated from single sinusoid. (b) Resulting 2D spectrum with 
energy peaks at correct wavenumbers. (c) Spectral cross section with vertical axis scaled 
by wavenumber. Note peak energy at correct wave number. 
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Figure 18. 2D autocorrelation plot for field created from a single sinusoid. Note the 
distortion caused by the angular binning process as angles come less perpendicular to the 
axis of the crests. 


2. Single Sine Wave Derived from Different Motion Data 


The second test case was another single sine wave with identical parameters as in 
Case 1. The object of this case was to determine if the output was independent of the 
velocity estimates used to translate the ULS time series into a spatial series. The identical 
field from Case 1 was generated but this time it was sampled by ice motion data for 
Mooring C from the previous year (August 2003-August 2004). Again, the spectra put 
the peak energy at the correct wave number, suggesting that the technique is insensitive 
to the details of the conversion of the time series to a spatial series with the buoy-derived 


velocities (Fig. 19). 
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Field from Single Sinusoid 
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Figure 19. (a) Field generated from sinusoid. (b) 2D spectrum from field sampled 
with different ice motion data. (c) Spectral cross section. 
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3. 2 Sine Waves of Different Wavelengths but Different Amps and 
Different Orientation 


This test case sampled a field made by combining two differing sinusoidal fields 
(Fig. 20). The first field had a wavelength of 200m, vertically oriented crests, and an 
amplitude of 5m. This corresponded to a wavenumber of k=.005 m” and log10(k)=-2.3. 
The second field had sinusoids of 500m, diagonally oriented crests, and amplitude of 


20m. This corresponded to a wavenumber of k=.002 m™ and log10(k)=-2.69. 


At a glance, it may appear that only the larger amplitude wave was represented in 
the 2D spectrum (Fig. 21(b)). However, cross sections through the 2 directions of 
interest show that both sinusoids were correctly resolved (Fig. 21(c)). Not only were 
peaks located in the correct wave numbers, but their relative sizes were consistent with 
their relative amplitudes. This reflects the larger variance in larger amplitude ice keels. 
This was one of the most important test cases because it showed that the process is robust 
enough to properly display a field created from sinusoids with differing orientation, 


wavelength, and amplitude. 
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Figure 20. Sinusoids #1 and #2 which combine to make up the composite field of two 
sinusoids with differing amplitude, wavelength, and orientation. 
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Composite Field from Two Sinusoids 
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Figure 21. (a) Composite field from two sinusoids. (b) 2D spectrum of field. (c) 


Spectral cross sections with peak energy and relative amplitude correctly represented. 
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4. Single Sinc Function 


This test case was to test how well the process worked on a non-sinusoidal field. 
A field was created with a sinc function repeating periodically with a wavelength of 
100m (Fig. 22(a)). To approximate the shape of a true ice keel, the amplitude was set at 
5m, it was 20m wide base. The peak amplitude of the spectra was at the correct 
wavenumber (Fig. 22(b)). However, harmonics were observed in the spectra (Fig. 22(c)). 
This is due to the nature of the 2D FFT process, which struggles to fit high wavenumber 
sinusoids into the observed profile. The harmonics were expected however, as the 
Fourier transform of the sinc function is the square wave. The implication is that the 


process can correctly represent non-sinusoids. 
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Single Sinc Function 
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Figure 22. (a) Field from single sinc function. (b) 2D spectrum from field. (c) 
Spectral cross section; note the decreasing harmonics at higher wavenumbers. 
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B. BEAUFORT GYRE DATA 


1. Mooring C: August 2004—August 2005 


Mooring C had the most dense buoy coverage and is used to illustrate the 
Beaufort Gyre data. The case described below is from August 2004- August 2005, the 
second year of deployment. This dataset had the best coverage with respect to the 2D 


autocorrelation function and therefore is considered the most robust. 


The 2D spectrum (presented as the square root of variance for display purposes) 
appears by eye to be directionally variable (Fig. 23a). Eight spectral cross sections taken 
in directions ranging from 000 to 157.5° in 22.5° intervals show that the spectra have a 
power law dependence on radius (Fig. 23(b) and (c)). Most of the variance of the 
underside morphology is contained in a wavenumber band between 0.001 and 0.01m" 
(100-1000m). Spectral levels ranged from 1.2 to 4.2 x 10‘ m*. There is no distinct, 


preferred wavelength. 
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Figure 23. Mooring C: 2004-2005. (a) 2D spectrum. (b) 8 energy conserving 
spectral cross sections (000 to 157.5° each 22.5°). (c) 8 spectral cross sections. 
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In order to determine whether or not the cross sections are statistically different 
from one another, the standard deviation needed to be determined. ‘This allowed 
confidence limits to be plotted against one another to check for overlapping regions, 
which would mean that there is not a difference. However, if the regions bounded by the 


confidence intervals were mutually exclusive, they are different and indicate anisotropy. 


Determination of the confidence intervals was not straightforward. Since the data 
was taken at one frequency, autocorrelated, subsampled, binned, meaned and gridded in 
autocorrelation space before passing through a 2D FFT, degrees of freedom are 
impossible to determine analytically. Therefore, the degrees of freedom were determined 
empirically by Monte Carlo simulation. A field was created which would ultimately 
approximate the characteristics of the observations. This was done by starting with a 
field of random numbers, taking the 2D FFT, multiplying by 1/k and then taking the 
inverse 2D FFT, to get back to a Cartesian field. This new field had the properties that 
when it was put through the process it yielded a 2D spectrum which was comparably red 


like the real data, but composed of random numbers. 


This random field was sampled using real buoy motion and put through the 
processing routine. A radial cross section was taken through the resulting spectrum. This 
entire process, from creating the random field through taking the cross section was 
repeated 30 times. For each wavenumber bin, the standard deviation of the cross sections 
was determined. A line was fit to approximate the standard deviation as a function of 


wavenumber. This represented the sigma value used for confidence limits. 


Two cross sections were taken in directions which appeared to be similar (Fig. 
24(a)). Even at 1 sigma, the spectra showed no regions of mutual exclusivity (Fig. 
24(b)). This suggested isotropy. At 2 sigma, the confidence intervals are even wider and 
the overlap even greater (Fig. 24(c)). The result was to conclude that the spectrum is 


isotropic in these two directions. 


Next, 2 cross sections were taken in directions which appeared to be significantly 
different (Fig. 25(a)). At 1 sigma, the spectra did show regions of mutual exclusivity, 


especially at the high wavenumbers, suggesting anisotropy (Fig. 25(b)). However, these 


4] 


regions of exclusivity disappeared at all except for a few points when the confidence 
intervals were set to 2 sigma, which suggests isotropy (Fig. 25(c)). The result was to 
conclude that there is ‘weak’ anisotropy in this spectrum. This is defined to mean that 


the results were isotropic except for low levels of statistical confidence. 
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Figure 24. (a) 2 radial slices through 2D spectrum. (b) | sigma confidence intervals 
for each slice. The red lines represent the + and — | sigma confidence intervals for the 
spectral slice taken along the red dotted line in (a) and the green lines are the confidence 
intervals for the slice taken along the green dotted line. (c) 2 sigma confidence intervals. 
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Figure 25. (a) 2 radial slices through 2D spectrum. (b) 1 sigma confidence limits for 
each slice. The red lines represent the + and — | sigma confidence intervals for the 
spectral slice taken along the red dotted line in (a) and the green lines are the confidence 
intervals for the slice taken along the green dotted line. (c) 2 sigma confidence intervals. 
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To ensure that the observed anisotropy was not a function of the binning process, 
a “red” field like the one used to generate the confidence intervals was sent through the 
process. Cross sections were taken in the 2 directions which showed the greatest 
amounts of anisotropy for the previous run. The result showed that the spectrum was 
isotropic, even at one sigma (Fig. 26). This shows that the observed isotropy was a 


characteristic of the real data and not the processing routine. 
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Figure 26. 1 sigma confidence intervals from the spectral cross sections from the 
synthetic “red” field taken in the same directions as Fig. 25. 


2. Summary for Other Data Sets 


There were directional differences observed in the spectra of all data sets at 1 
sigma. These suggest anisotropy. However, this was never true at 2 sigma, where the 
spectra appeared isotropic. Furthermore, although the isotropic regions observed usually 
occurred at high wavenumbers, they were also sometimes observed throughout the rest of 
the spectra. As for the data from Mooring C from 2004-2005, most of the variance of the 
other datasets is contained in a wavenumber band between 0.001 and 0.01m' (100- 
1000m). Again, there are no distinct, preferred wavelengths. Spectral levels ranged from 
0.2 - 1.6 x 10 m? for Mooring A from 2004-2005 to 0.3 - 3.3 x 10* m? for Mooring C 
from 2003-2004. 
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In summary, the spectra are weakly anisotropic for the wavenumbers resolved in 


this study. 
C. SPATIAL AND TEMPORAL COMPARISONS 


The weak anisotropy observed was not consistent from one mooring to another. 
The 2D spectra are displayed in Figs. 27 and 28. Spatial and temporal comparisons failed 
to determine a consistent directional preference. However, most of the variance was 
contained in the .001 to .01m! wavenumber band (100-1000m) for all datasets. Spectral 


levels varied significantly between the datasets. 
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Figure 27. 2D spectra of under ice roughness from August 2003-August 2004. (a) 
Mooring A. (b) Mooring B. (c) Mooring C. 
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Figure 28. 
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Mooring A. (b) Mooring B. (c) Mooring C. 
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D. COMPARISONS TO OTHER RESEARCHERS 


This weak anisotropy does not agree with the findings of Davis and Wadhams 
(1995). While working from a submarine, they gathered concurrent upward looking 
echosounder and towed sidescan sonar which allowed for determination of both ice draft 
and ice keel orientation. Their findings showed that there was anisotropy in the 
alignment of the large ice roughness elements, the ice keels. Their data showed lineation 
which fell roughly into two preferred alignments. When considered in the context of the 
area in which the data were gathered, the directions are roughly along and perpendicular 
to the Greenland coast (Fig. 29). The observed ice keels could be the result of ridging 
due to compression caused by downstream boundary and shearing caused by motion of 


Transpolar Drift Stream. 
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Figure 29. Ridge orientation vectors (Davis and Wadhams, 1995). 


Comparisons to the findings of Davis and Wadhams and this project should be 
made carefully. This project used data from the Beaufort Gyre, well away from land 
boundaries. Furthermore, the data was collected from moorings in deep water (greater 
than 3500m) so there was no bottom interaction or grounding to cause a preferred angle 


of orientation for roughness elements. 
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Hibler (1972) first suggested the process upon which this project was based. 
Using his version of data processing, he found anisotropy to suggest a preferred lineation 
of pressure ridge sails in the Beaufort. However, this was based on one case using three 
intersecting laser profiles about 2400m long collected by aircraft. In contrast, this project 
used several hundred profiles in three locations in the Beaufort over the course of two 


years. Again, comparisons must be made carefully. 


E. TURBULENCE IMPLICATIONS 


The implications of the weak anisotropy that was observed is that directionality is 
important at very small scales. In turbulence simulations, where roughness is significant, 
a ‘spectral’ approach is most appropriate to approximate a near-isotropic field. In 
contrast, most conceptual models of flow over roughness have been developed with 
regular arrays of ridges and keels with a single orientation and fixed spacing (Arya, 


1975). 


Density profiles obtained at APLIS-2007 in the Beaufort Sea showed pycnoclines 
which were significantly shallower than those historically observed (Fig. 30). This 
shoaling of the thermocline is possibly explained by the recent accelerating rate of ice 
melt. Pycnoclines at shallow depths can produce the “dead water” effect if an ice keel is 
dragged through them. With any anisotropy in the ice roughness, especially for deep ice 
keels whose depth is comparable to or exceeding the pycnocline depth, there will be 
significant effects on the momentum budget. Energy will be lost at variable rates 


depending on the direction that the ice flow is drifting. 
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Figure 30. Density profiles from APLIS 2007. Note the shallow pycnocline; 
historically the pycnocline has been observed near 35m. 


F. POSSIBLE MITIGATING EFFECT 


An observation made recently showed that the ice floe used for APLIS-2007 
rotated significantly, changing heading by 11 degrees in 7 days (Fig. 31) (Hutchings, 
2007). This has the potential to mitigate the effects of weak directional dependence that 
were observed in this study, making the time averaged roughness isotropic. The APLIS- 
2007 observation is considered abnormal when compared to previous ice camp rotation 
measurements. However, as the Arctic pack ice becomes thinner and more mobile, more 
rotation is possible. This would make conditions more isotropic, and should be 


investigated further to determine if it is a new trend. 
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Figure 31. Plot of APLIS 2007 ice floe rotation (Hutchings, 2007). 
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Vv. CONCLUSIONS 


The roughness of the ice underside plays a role in turbulent ocean heat, salt and 
momentum fluxes that determine the delicate balance between surface forcing and the 
ocean interior that sustains or reduces the perennial ice cover. It has climate, 


engineering, and naval implications as well. 


This project took sparsely sampled ice draft data and processed it in such a way 
that the 2D wavenumbers were conserved in order to produce a realistic 2D spectrum of 
under-ice roughness. The process was validated using several test cases. It showed weak 
anisotropy of under-ice roughness at scales between 10m and 2000m for pack ice in the 
Beaufort Gyre. Error estimates were determined by Monte Carlo simulations. For each 
dataset, most of the variance of the underside morphology was contained in a 
wavenumber band between 0.001 and 0.01m' (100-1000m). There were no distinct, 
preferred wavelengths. Spectral levels ranged from 1.2 - 4.2 x 10° m? for Mooring C 


(2004-2005) to 0.2 — 1.6 x 10* m? for Mooring A (2004-2005). 


Future studies will benefit from the new mooring configurations in the Beaufort 
Gyre Exploration Project. The new moorings have colocated ULS and ADCPs which 
when set in bottom tracking mode give the true local ice motion at the mooring. There no 
longer needs to be interpolation between buoys or assumptions of continuous ice motion. 
This increased amount of data will allow tighter limits in the processing routine and 


ultimately a higher resolution 2D roughness spectrum. 
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